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We investigate the gravitational evolution of dark matter halos made up of a massless bosonic 
field. The coupled Einstein-Klein-Gordon equations are solved numerically, showing that such a 
boson halo is stable and can be formed under a large class of initial conditions. We also present an 
analytical proof that such objects are stable in the Newtonian limit. In the context of boson stars 
made of massive scalar fields, we introduce new solutions with an oscillatory scalar field, similar to 
boson halos. We find that these solutions are unstable. 
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I. INTRODUCTION 

X) ■ 

^ I Boson stars are self-gravitating configurations which are described mathematically by the Einstein-Klein-Gordon 
equations . In order to find such globally regular solutions it is necessary to have a conserved particle number which 
■ is fulfilled if one uses a complex scalar field and a potential with global U{1) symmetry. Investigations so far have 
shown that stable and unstable configurations exist so that study to a greater extent makes sense . A boson star 
can have values of the mass over large orders of magnitudes; it depends on the mass of the scalar field and possible 
^ [ self- interactions ||^ . A recent paper described the effects if the scalar field would interact only gravitationally so that 
the star is transparently similar as a dark matter halo of a galaxy |^]. Baryonic matter could then move within the 
gravitational potential of the star and if perhaps a accretion disk can be built the radiation would suffer a gravitational 
. redshift. Redshift of an iron-K line is actually observed within the center of a Seyfert galaxy Q. 
' Boson star solutions are characterized by an exponential decrease of the scalar field for which a mass term in the 
^ potential is responsible. A zero temperature solution is usually considered in which all scalar particles are within the 
' same ground state forming a Bose-Einstein condensate. The size of this general-relativistic object is given by the 
■S^^ ' Compton wavelength. If one 'switches off' the mass of the scalar field, the Compton wavelength is infinite. Then 
O , the massless scalar particles are within a coherent state, a kind of 'boson star', which in principle has an infinite 
^ '' range. Investigations of these solutions which we call now boson halos, have shown that the mass of these 

?H I objects increases linearly, a behavior suggested by HI observations for spiral and dwarf galaxies. It was shown that 
the rotation curves for galaxies could be very well fitted by this kind of matter. In this paper, we present an analytical 
J> , proof showing the stability of the Newtonian solutions. 

In this paper we also examine numerically how these objects could have formed and show that they are stable. For 
rN this purpose, we use the process of gravitational cooling as described in ||l^. We add a small radial perturbation to 
a solution and show that the evolved solution does not disperse or form a black hole, hence it is stable. Furthermore, 
we present here for the first time solutions with an oscillating scalar field for the massive case; we find that they are 
unstable. 

In section II we present the mathematical foundations of the problem of bosonic objects made from massless scalar 
fields. We set up equilibrium equations and approximate solutions. We then set up the evolution equations which we 
use in our code. In section III we discuss our numerical studies on the formation of scalar halos followed in section 
IV by an analytical proof of stability of these objects. In section V we discuss our numerical studies of stability of 
bosonic halos. In section VI we discuss new oscillatory (in space) solutions for massive scalar fields. Finally we have 
a discussion of our results and their implications in section VII. 
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II. MATHEMATICAL FOUNDATIONS OF THE PROBLEM 



In this Section we start with the Lagrangian for the system and set up the equiUbrium equations and their solutions 
as well as the evolution equations with which we evolve the equilibrium solutions. 



A. Einstein-scalar-fleld equations 



The Lagrange density of a complex self-gravitating scalar field reads 



2 V I y 



i? + 5''^(a^$*)(5,$)-i7(|$|2) 



(1) 



where R is the curvature scalar, k = 8nG, G the gravitation constant {h = c = 1), g the determinant of the metric 
g^^, <I> the complex scalar field, and U the potential depending on |$|^ so that a global U{1) symmetry is conserved. 
Then we find the coupled system 



where 



is the energy-momentum tensor and 



a = d„ 



9 \9"'d. /v17 



the generally covariant d'Alembertian. 

For spherically symmetric solutions we use the following line element 

ds^ = N^{r,t)de - g^{r,t)dr^ - r^{di&^ + sin^ ^dip^) 

and for the scalar field the ansatz 

$(r,t) =P(r,t)e-''"* , 

where u) is the frequency. 
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B. Equilibrium Equations 



The non-vanishing components of the energy-momentum tensor are 



-T2'' = -T3^=P± = 



^2p2 p/2 



3/2 
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^2p2 p,2 
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where ' = d/dr. 

The equilibrium configurations (those for which the metrics are static) for a system of massless scalar fields are 
derived from the following equations 
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(11) 
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Here we use the dimensionless variables r = ojr, t = ujt and a = \/AttG (j). Regularity at the center imphes g{r = 0) = 1. 
The equihbrium configurations are characterized by saddle points in the density p, which might have some significance 
in stability issues. These systems are characterized by a two parameter family of solutions. Firstly, we can vary a{0), 
the central density, and secondly, for each a{0) we can vary A^(0). Although we do not have asymptotic flatness (the 
mass M rises with r) the mass profile for a given value of r has a very interesting feature. Figure la shows the mass 
as a function of central density for a given value of r = i? and N(0). The mass throughout the paper is calculated 
from the Schwarzschild metric 



M(r) = 



1 - 



g- 



(15) 



where g^ is the radial metric. 

The mass profile is very similar to that of the massive scalar field case although there the mass is not subject to 
these variations with r. This is also reminiscent of neutron star profiles. The increase in mass with increase in central 
density is followed by a decrease in mass with further increase in central density. For the same value of N{0) and 
larger values of r the mass increases but the profile remains the same. Figure lb shows mass versus radius for different 
central densities and a given value of N(0). The profile is independent of the value of the radius for large values of r. 

The Noether theorem associates with each symmetry a locally conserved "charge". The Lagrangian density (|l|) 
is invariant under a global phase transformation $ — s- $e^"'. The local conservation law of the associated Noether 
current density reads 



df.f = 



(16) 



If one integrates the time component j° over the whole space we find the particle number 



Np = Attuj J j^r^P^ dr 





(17) 



C. Approximate solution 



The Newtonian and the general relativistic solutions of system (^T])-(^4|) for [/ = were recently discussed [p|-pT|. 
The Newtonian ones (initial value cr(0) smaller than about 10~^) can be used to describe dark matter halos of spiral 
and dwarf galaxies. The general relativistic ones (initial value a(0) > 10~^) reveal very high redshift values. As 
was shown in [p|-pl|, for the Newtonian solutions there exists an approximately analytic formula for the scalar field 
{N^g^l) 



a{r) = A 



sin(r) 



(18) 



where A is a constant. The energy density reads 



p(r) 



1 



sin(2r) sin^(r) 



(19) 



and the mass defined as usual as M{r) — piC)C'^d( yields 
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M(r) = 



cos(2r) — 1 



(20) 



The energy density shows a decreasing behavior with saddle points in between (cf. Fig. 3a). The Newtonian sohition 
for the particle number is 
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Mr) = -T 



2r - sin(2r) 



(21) 



the quantitative behavior of which is also revealed by numerical calculation. 



D. Evolution Equations 



For accurately numerical evolution the following set of variables are chosen 



TTl 



1 9^1 
a dt 



T^2 



1 04^2 

a dt 



where 



N 
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(22) 



(23) 



and the subscripts on ipi denote the real and imaginary parts of the scalar field multiplied by r. 

In terms of these variables and the dimensionless ones in the previous Section the evolution equations are as follows: 
The radial metric function g evolves according to 



g = N[TTia[ + 7r2cr2) . 
The polar slicing equation, which is integrated on each time slice, is given by 

N 



N' = 



+ r[{a[r + {a'2r]+^ 



The Klein-Gordon equation for the scalar field can be written as 



TTi — a'4''i + aip'/ — Tpi 



gN+- 
r 



z = l,2, 



(24) 



(25) 



(26) 



ipi = a-Ki, z = 1, 2 . 
The Hamiltonian constraint equation is given by 



1 TT 



t'2 I ^'2 



J-2g2 



^2g2 







(27) 



(28) 



After we introduce a perturbation in the field or when we have an arbitrary initial field configuration we reintegrate 
on the initial time slice to get new metric components B. 



III. FORMATION OF SCALAR OBJECTS 



In this Section we discuss the possibility of formation of massless scalar field configurations from non-equilibrium 
data. Earlier studies both numerical and analytical have confirmed that self-gravitating objects with massless scalar 
fields cannot be compact |12 1^. Thus taking a Gaussian type localized distribution of scalar matter and evolving the 



system using the evolution equations described above, results in dissipation of the scalar matter without forming any 
self-gravitating object. Even sinusoidal functions that were damped by exponential decays met with the same fate. 
These were functions like exp(— r) sin(r)/r and exp(— r) cos(r). In Figure 2, a plot of the density versus the radius for 
different times is shown for a = 0.001 cos(r) exp(— r). The star dissipates very quickly. 
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The initial configurations for the scalar field that yield stable configurations were those characterized by a l/r times 
a sinusoidal dependence at large r as well as a saddle point structure of the energy density p. These were typically 
functions like 



cos(r) 



cos(r) 



1 — exp 



and 



(29) 



sin(r) 1 + - log 



1 



(30) 



In Figure 3 we show the density and radial metric evolutions for a field of the form 0.003 cos(r)(l — exp(— l/r)). This 
settles into a self-gravitating object after some time. Figure 3a shows the density as a function of r. The central 
density p(0) increases during the evolution from its initial value at t = 0. Figure 36 reveals the radial metric as it 
evolves in time as a function of radius. This too displays the settling to a stable configuration. In Figure 3c we show 
the mass loss for the system as it finds itself a configuration. The mass at fixed radial values for these times is shown 
in Table 1. The amount of radiation for this system is relatively small and decreases in time. 

On the other hand, functions like l/(r + 1) , cos(r)/(l + r)^ and sin(r^ °^)/r^ *'^) failed to settle to a bound state 
and just dissipated away. These functions did not have saddle points in p. 

Functions like a = sin(r) / log(l + r) for which p has saddle points (it also has a maximum near the origin), did not 
dissipate but failed to settle in a long numerical evolution. This might mean they would form the halo structures in 
a very long time for which it would not be computationally feasible to evolve the system. The same was observed for 
sin(r)/(r' log(l + r)) or a numerical six node boson star configuration [uj < m). 

All the above configurations were fairly Newtonian for which the exact solution was close to sin(r)/r. 



IV. ANALYTICAL STABILITY PROOF 



We investigate here the stability of the Newtonian solutions against small radial perturbations; for stability inves- 
tigations concerned with the boson star model and using a perturbative method, see Ref. [p^-p^. Therefore, we can 
neglect the perturbations for the spacetime and make the following ansatz for the linear scalar field perturbations SP: 

$(r, t) := P(r)e-'"* -f (5F(r)e*'="* , (31) 

where kn are the frequencies of the normal modes. The pulsation equation, an eigenvalue equation, is then: 

SP" + 26P'/x + kl5P = . (32) 

Together with the boundary conditions 

5P'{Q) = 5P{R) = (33) 

for regularity of 5P at the origin and at the radius R of the boson halo, this is a Sturm-Liouville eigenvalue problem. 
As it is well-known from mathematical theory, we have a series of real eigenvalues with a minimal one: 

kl<kl< ... . (34) 



Eigenfunctions of this particular differential equation (|3^) exist only if the eigenvalues k'^ are positive. This means 
that all modes are stable. The eigenfunctions are 



where the eigenvalues are 



5P , (35) 



kn = n^ ; (36) 



n is a natural number and R the radius of the solution. 
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V. NUMERICAL STABILITY ISSUES 



In order to study the stability of these systems numerically we perturbed exact configurations and observed whether 
the system settled to new configurations. In some sense taking a function like As'm{r)/r which is quite close to the 
exact solution for a Newtonian system or A cos{r) / {1 + r) can be regarded as a perturbation on the exact solution and 
the eventual settling to a new configuration can be regarded as indicative of stability. The settling of 0.001 cos(r) / 
to a stable configuration is shown in Figure 4. The density at the center increases in this case from what it starts at. 
The radial metric evolution case is shown in Figure 4b. On the other hand when A = 1 the system is non-Newtonian 
and in a long evolution failed to settle down although it did not disperse. This is not surprising since the sin(r)/r 
form is close to the solution only in the Newtonian case. 

In Figure 5 we show a perturbed Newtonian configuration as it settles to a new Newtonian configuration. The 
perturbation that is used in this case mimics an annihilation of particles. A Gaussian bump of field is removed from 
a part of the star near the origin. Figure 5a is a plot of the unperturbed versus the perturbed density at t = 0. The 
perturbed configuration is evolved and settles to a new configuration. The scalar radiation moves out as shown in 
Figure 5b. In Figure 5c the density profile is shown after the system settles down. The system is very clearly in a 
new stable configuration. In Figure 5d the mass is plotted as a function of radius for various times. Again one can 
see that the mass loss is decreasing by the end of the run showing that the system is settling to a new configuration. 
The mass as a function of time is presented in Table 2 for different radii. 

We have so far been successful only in our Newtonian evolutions. The reason for this is that denser configurations 
need much better resolution for evolution. This is coupled with the difficulty that we still need the boundary to be 
very far away so that the density has significantly fallen off. We are working on improving our boundary conditions. 
So far, we are using either an outgoing wave condition or exact boundary conditions where the latter one simulates a 
vacuum energy. Further investigation is needed before we can decide whether the non- Newtonian configurations are 
inherently unstable. 

VI. BOSON STAR OSCILLATORS 

A boson star consists of scalar massive particles, hence in the simplest model one has a potential U — r7i^|<I>p. 
Exponentially decreasing solutions exist for special eigenvalues ui < m, so that the star has a finite mass. In the case 
of w > m, oscillating scalar field solutions can be found for all values of u. The energy density reveals minima and 
maxima as opposed to the saddle point structure seen for the massless case. Figure 6a shows a comparison of the 
density profile for the two cases. 

In figure 6b we show the mass profile and the particle number as functions of central density. The binding energy is 
always positive making the system unstable against a collective transformation in which it disperses into free particles 
but the system is still stable against any one-at-a-time removal of a particle [T^ . The complete dispersion is a 
numerical confirmation of the result by Zeldovich, cf. [Q. 

VII. DISCUSSION 

We have shown that massless complex scalar fields are able to settle down to a stable configuration. It seems 
that the formation process needs a special form for the energy density. The appearance of saddle points within the 
decreasing density supports the 'birth' of the boson halos. In contrast, extremal values for the density cannot be 
compensated and lead to the destruction of the initial configuration. This result of our numerical investigation reveals 
that the formation of boson halos could be difficult. But as our attempt by cutting off a part of a Newtonian solution 
has shown, if there exists such a structure, then a new configuration forms and settles down to a stable dark matter 
halo. One can understand this as Bose-Einstein condensation where a small initial configuration forms and particles 
in the surrounding falls into the most favorable state given by the condensation. In this way, the boson halo forms 
and grows up to the point where it is stopped by the outer boundary, the vacumm energy density. Our choice of the 
numerical boundary was necessary so that it simulates a hydrostatic equilibrium with this vacuum energy as it was 
assumed in Because the boundary can be placed at every radius the size of these dark matter halo depends 

on the value of the cosmological constant. 

During the formation process of the boson halo, massless scalar particles change into the Bose-Einstein condensate 
and loose energy through this process; only particles with at least this energy can participate. We infer that the 
condensed state is the energetically more preferable state in comparison with the free state. The calculation of the 
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binding energy, the comparison of mass of the bounded with the free particle at infinity, is not advised due to the 
arbitrariness of the energy of a massless particle at infinity. 

Investigations with real scalar fields |13| or with different boundary conditions ^] cannot lead to singularity free 
solutions; cf. Liddle and Madsen in 

The second class of solutions presented here is the oscillatory one for a massive complex scalar field. Our result of 
a positive binding energy confirms the expectations for systems where the energy of the field is higher than the rest 
mass energy. Our numerical evolution also reveals the instability by dispersion of the unperturbed exact solution. 
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FIG. 1. (a) Left: The mass profile of equilibrium configurations of massless scalar fields, are interestingly similar to the 
profiles of massive scalar fields and neutron stars. The mass for a given radius increases to a maximum with central density 
before it decreases with further increase in central density. However by calculating the particle number, we see no division into 
stable and unstable configurations at the peak. The particle number depends very sensitive on the initial value of the lapse 
function A'', so that a clear cusp structure in the mass-particle-number diagram found for boson stars cannot be derived here 
(cf. (b) Middle: The profiles are largely independent of radius although the mass itself increases with r. (c) Right: If 

one cuts the boson halos at the same energy density, then the mass profile reveals the same oscillating behavior as in case (a). 
The energy density p is measured in up' jn. 




FIG. 2. The inability of massless scalar field configurations to form compact self-gravitating objects as discussed by 
Christodoulou and others is verified numerically. An initial field configuration of the form a = 0.001 cos(r-) exp(— r) is 

seen dispersing in the plot. The dispersion takes place very quickly. 
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FIG. 3. (a) Left: The formation of a self-gravitating massless scalar field object is shown here for an initial field configuration 

of the form 0.003cos(r)(l — oxp( — 1/r)). The density in this case increases from its initial value as the configuration settles 
down, (b) Middle: The radial metric for this configuration is shown evolving to a stable final configuration, (c) Right: The 
mass is plotted as a function of radius at different times. As the configuration evolves it looses less and less mass as it settles 
down. 




FIG. 4. (a) Loft: A Newtonian initial field configuration of the form 0.001 cos(r)/(r + 1) can be regarded in some sense to 
be a perturbation of a Newtonian configuration which is proportional to sin(r)/r. A plot of the density evolution is shown 
as it settles to a final stable state. Configurations of this form with non-Newtonian amplitudes do not settle down. Since 
the non-Newtonian configurations do not have this sinusoidal dependence this is not surprising, (b) Right: The radial metric 
evolution for the above configuration shows its settling in time. 
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FIG. 5. (a) Above left: An exact Newtonian configuration is perturbed and evolved. The initial unperturbed and perturbed 
densities are shown. A Gaussian amount of scalar field is removed fi'om the equilibrium configuration and the constraint 
equations axe reintegrated to give new metrics before beginning evolution, (b) Above right: The outgoing scalar radiation can 
be seen in a plot of the density at early times, (c) Bottom left: The density profile is shown after it starts to settle down. 
Between times of 300 and 1500 very little has changed in the density profile, (d) Bottom right: The mass is plotted as a 
function of radius for different times. In order to enhance the features only a small part of the radial region is shown. The 
amount of mass loss is clearly lessening in time. 
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FIG. 6. (a) Left: The density profiles of massive scalar field configurations with cu > m is compared to that of a masslcss 
field. The former has maxima and minima while the latter has saddle points. This structure may be related to stability issues, 
(b) Right: The mass and particle number are plotted against the central density for boson star oscillators. We find in all cases 
that the mass is always greater than the particle number explaining their instability against a collective dispersion. 



TABLE L Mass values at different radii and different times for the evolution in Fig. 3. One recognizes that the change in 
the mass at fixed radius decreases, hence it proves the settling to a stable configuration. 



Table 1. 





t = 


t = 500 


t = 700 


t = 900 


r = 800 


0.03589375 


0.0350625 


0.034625 


0.034234375 


r = 1200 


0.0575 


0.0563125 


0.055625 


0.0550625 


r = 1600 


0.0801875 


0.078625 


0.0776875 


0.0769375 



TABLE II. Mass values at different radii and different times for the evolution in Fig. 5. One recognizes that the change in 
the mass at fixed radius decreases, hence it proves the settling to a stable configuration. 



Table 2. 





t = 


t = 1500 


t = 2500 


t = 3400 


r = 500 


0.002847 


0.00278 


0.0026875 


0.00268375 


r = 1000 


0.00628625 


0.00619375 


0.00608105 


0.00608335 
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